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1 Introduction 

Fiber bundles with statistically distributed thresholds for breakdown of indi- 
vidual fibers are interesting models of the statics and dynamics of failures in 
materials under stress. They can be analyzed to an extent that is not possi- 
ble for more complex materials. During the rupture process in a fiber bundle 
avalanches, in which several fibers fail simultaneously, occur. We study by 
analytic and numerical methods the statistics of such avalanches, and the 
breakdown process for several models of fiber bundles. The models differ pri- 
marily in the way the extra stress caused by a fiber failure is redistributed 
among the surviving fibers. 

When a rupture occurs somewhere in an elastic medium, the stress else- 
where is increased. This may in turn trigger further ruptures, which can cas- 
cade to a final complete breakdown of the material. To describe or model such 
breakdown processes in detail for a real material is difficult, due to the com- 
plex interplay of failures and stress redistributions. Few analytic results are 
available, so computer simulations is the main tool (See Refs. [1], [2] and [3] for 
reviews). Fiber bundle models, on the other hand, are characterized by simple 
geometry and clear-cut rules for how the stress caused by a failed element is 
redistributed on the intact fibers. The attraction and interest of these models 
lies in the possibility of obtaining exact results, thereby providing inspiration 
and reference systems for studies of more complicated materials. 

In this review we survey theoretical and numerical results for several mod- 
els of bundles of N elastic and parallel fibers, clamped at both ends, with 
statistically distributed thresholds for breakdown of individual fibers (Fig. 1). 
The individual thresholds Xi are assumed to be independent random variables 
with the same cumulative distribution function P{x) and a corresponding den- 
sity function p{x): 



Prob(a;i < a;) = P{x) = / p{u) du. 

Jo 



(1) 
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Fig. 1. A fiber bundle of N parallel fibers clamped at both ends. The externally 
applied force is F. 

Whenever a fiber experiences a force equal to or greater than its strength 
threshold Xi , it breaks immediately and docs not contribute to the strength of 
the bundle thereafter. The maximal load the bundle can resist before complete 
breakdown of the whole bundle is called the critical load. The models differ 
in the probability distribution of the thresholds. Two popular examples of 
threshold distributions are the uniform distribution 

pw={t' If (^) 

and the WeibuU distribution 

P(x) = 1 - eyip{-{x/xr)''). (3) 

Here .t > 0, reference threshold and the dimensionless number k is the 

WeibuU index (Fig. 2). 
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Fig. 2. The uniform distribution (A) and Weibull distributions (B) with fe = 5 (solid 
line) and A; = 10 (dotted line). 



Much more fundamental, however, is the way the models differ in the 

mechanism for how the extra stress caused by a fiber failure is redistributed 
among the unbroken fibers. The simplest models are the equal-load-sharing 
models, in which the load previously carried by a failed fiber is shared equally 
by all the remaining intact bonds in the system. That some exact results 
could be extracted for this model was demonstrated by Daniels [4] in a classic 
work some sixty years ago. Local-load-sharing models, on the other hand, are 
relevant for materials in which the load originally carried by a failed fiber is 
shared by the surviving fibers in the immediate vicinity of the ruptured fiber. 

The main property of the fiber bundle breakdown process to be studied in 
the present review is the distribution of the sizes of the burst avalanches. The 
burst distribution D{A) is defined as the expected number of bursts in which 
A fibers break simultaneously when the bundle is stretched until complete 
breakdown. For the equal-load-distribution models that we consider in Sec. 2 
Hemmer and Hansen [5] showed that the generic result is a power law, 

lini ^ oc (4) 

for large A, with ^ = 5/2. In Sec. 2.2 we will show, however, that for some 
rather unusual threshold distributions the power law (4) is not obeyed. More 
importantly, we show in Sec. 2.4 that when the whole bundle at the outset is 
close to being critical, the exponent ^ crosses over to the value ^ = 3/2. In 
Sec. 2.5 we pay particular attention to the rupture process at criticality, i.e., 
just before the whole bundle breaks down. 
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The average strength of the bundle for a given load can be viewed as 
the result of a sequential process. In the first step all fibers that cannot with- 
stand the applied load break. Then the stress is redistributed on the surviving 
fibers, which compels further fibers to burst. This starts an iterative process 
that continues until equilibrium is reached, or all fibers fail. When equilib- 
rium exists, it characterizes the average strength of the bundle for the given 
load. This recursive dynamics can be viewed as a fixed-point problem, with 
interesting properties when the critical load is approached. We review such 
recursive dynamics in Sec. 2.6. 

For other stress redistribution principles than equal-load-sharing, the 
avalanche distributions arc different from the power law (4). In Sec. 3 we 
study examples of such systems. Special cases are local-stress-distribution 
models in which the surviving nearest neighbors to a failed fiber share all the 
extra stress, and a model in which the fibers are anchored to an elastic clamp. 

2 Equal-load-shciring fiber bundles 

This is the fiber-bundle model with the longest history. It was used by Pierce, 
in the context of testing the strength of cotton yarn [6] . The basic assumptions 
are that the fibers obey Hookean elasticity right up to the breaking point, and 
that the load distributes itself evenly among the surviving fibers. The model 
with this democratic load redistribution is similar to mean-field models in 
statistical physics. 

25 I ■ . ■ ' ■ 1 
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Fig. 3. F(x) vs. X curve. Avalanches are shown as horizontal lines. 



At a force x per surviving fiber, the total force on the bundle is 

F{x) = Nx[l - (pix)] , (5) 



Rupture processes in fiber bundle models 



5 



where (t>{x) is the fraction of failed fibers. In Fig. 3 we show an example of 
a F vs. X. We have in mind an experiment in which the force F, our control 
parameter (Fig. 1), is steadily increasing. This implies that not all parts of 
the F{x) curve are physically realized. The experimentally relevant function 
is 

Fp^a;) = LMF , (6) 

the least monotonic function not less than F{x). A horizontal part of Fph{x) 
corresponds to an avalanche, the size of which is characterized by the number 
of maxima of F{x) within the corresponding range of x (Fig. 3). 

It is the fluctuations in F(x) that create avalanches. For a large sample 
the fluctuations will be small deviations from the average macroscopic char- 
acteristics {F). This average total force is given by 

{F){x)=Nx[l-P{x)\. (7) 

Let us for the moment assume that {F) {x) has a single maximum. The maxi- 
mum corresponds to the value x = Xc for which d{F) jdx vanishes. This gives 

1 - P{Xc) - Xc pixc) = 0. (8) 

The threshold Xc corresponding to the maximum of F is denoted the critical 
threshold. Because of fluctuations, however, the maximum value of the force 
may actually occur at a slightly different value of x. 



2.1 Burst distribution: The generic case. 
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Fig. 4. The burst distribution D{A)/N for the uniform distribution (A) and the 
Weibull distribution with index 5 (B). The dotted lines represent the power law 
with exponent —5/2. Both figures are based on 20000 samples of bundles each with 
N = 10^ fibers. 



6 



Per C. Hemmer, Alex Hansen, and Srutarshi Pradhan 



That the rupture process produces a power-law decay of the burst distri- 
bution D{A) is seen at once by simulation experiments. Fig. 4 shows results 
for the uniform threshold distribution (2) and the Weibull distribution (3) 
with index k = 5. 

In order to derive analytically the burst distribution, let us start by con- 
sidering a small threshold interval {x, x + dx) in a range where the average 
force {F){x) increases with x. For a large number N of fibers the expected 
number of surviving fibers is A'^[l — P{x)]. And the threshold values in the 
interval, of which there are Np{x)dx, will be Poisson distributed. When N is 
arbitrarily large, the burst sizes can be arbitrarily large in any finite interval 
of X. 

Assume that an infinitesimal increase in the external force results in a 
break of a fiber with threshold x. Then the load that this fiber suffered, will 
be redistributed on the N[l — P{x)] remaining fibers; thus they experience a 
load increase 

= N[i-Pix)y 

The average number of fibers that break as a result of this load increase is 

a = a{x)= Np{x) ■ 5x = ^^j^ ■ (10) 

For a burst of size A the increase in load per fiber will be a factor A 
larger than the quantity (9), and an average number a{x)A will break. The 
probability that precisely A — 1 fibers break as a consequence of the first 
failure is given by a Poisson distribution with this average, i.e. it equals 

e-- (11) 

{A-iy. • ^ ' 

This is not sufficient, however. We must ensure that the thresholds for these 
A — 1 fibers are not so high that the avalanche stops before reaching size A. 
This requires that at least n of the thresholds are in the interval (x, x + nSx), 
for 1 < n < Zi — 1. In other words, if we consider the A intervals (x, x + 6x), 
{x + 6x, X + 2Sx), . . . , {x + (A — 1)S, x + A5x), we must find at most n — 1 
thresholds in the n last intervals. There is the same a priori probability to find 
a threshold in any interval. The solution to this combinatorial problem is given 
in Ref. [7] . The resulting probability to find all intermediate thresholds weak 
enough equals \/A. Combining this with (11), we have for the probability 
(j){A, x) that the breaking of the first fiber results in a burst of size A: 

<P{A, x) = a(a;)^-ie-«(-)^. (12) 

This gives the probability of a burst of size Z\, as a consequence of a fiber 
burst due to an infinitesimal increase in the external load. However, we still 
have to ensure that the burst actually starts with the fiber in question and 
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is not part of a larger avalanche starting with another, weaker, fiber. Let us 
determine the probability Pb{x) that this initial condition is fulfilled. 

For that purpose consider the d—1 fibers with the largest thresholds below 
X. If there is no strength threshold in the interval {x — Sx,x), at most one 
threshold value in the interval (x — 2Sx, x), ... , at most d — I values in the 
interval {x — dSx,x), then the fiber bundle can not at any of these previous 
a;-values withstand the external load that forces the fiber with threshold x 
to break. The probability that there are precisely h fiber thresholds in the 
interval {x — Sx d, x) equals 

hi 

Dividing the interval into d subintervals each of length Sx, the probability Ph,d 
that these conditions are fulfilled is exactly given by ph,d = 1 — h/d (See Ref. 
[7]). Summing over the possible values of h, we obtain the probability that 
the avalanche can not have started with the failure of a fiber with any of the 
d nearest-neighbor threshold values below x: 

^^(-M) = E ^ e-^'i^--,) = (l-)e-''E ^ + ^ e-"^- (13) 

h=0 h=0 

Finally we take the limit d oo, for which the last term vanishes. For a > 1 
the sum must vanish since the left-hand side of (13) is non-negative, while the 
factor (1 — a) is negative. For a < 1, on the other hand, we find 

Pbix) = Km Pb{x\d) = 1 - a, (14) 

d — ^oo 

where a = a{x). The physical explanation of the different behavior for a > 1 
and a < 1 is straightforward: The maximum of the total force on the bundle 
occurs at Xc for which a{xc) = 1, see Eqs. (8) and (10), so that a{x) > 1 cor- 
responds to X values almost certainly involved in the final catastrophic burst. 
The region of interest for us is therefore when a{.x) < 1, where avalanches 
on a microscopic scale occur. This is accordance with what we found in the 
beginning of this section, viz. that the burst of a fiber with threshold x leads 
immediately to a average number a{x) of additional failures. 

Summing up, we obtain the probability that the fiber with threshold x is 
the first fiber in an avalanche of size A as the product 

^x) = cf>{A,x)Pb{x) = a(x)^-ie-''(-)^[l - a{x)], (15) 

where a{x) is given by Eq. (10), 

X p{x) 



a{x) 



1 - P{x) 



Since the number of fibers with threshold values in (a;, x+6x) is N p{x) dx, 
the burst distribution is given by 
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^ = ^ £ ^(^)^^(^) = ^ £ a{x)^-^e-<^^^ [1 - a{x)] p{x) dx. 

(16) 

For large A the maximum contribution to the integral comes from the 

neighborhood of the upper integration limit, since a(.x) e^"*^^) is maximal for 
a{x) = 1, i.e. for x = Xc- Expansion around the saddle point, using 

a^e-"^ = exp [A{-1 - \(l - af + 0{l - af)] , (17) 

as well as a{x) ~ 1 + a'{xc){x — Xc), produces 

^ = ^^:^a\x.)£p{x.) e-'(-)^(--)^^/^x-a.e) dx. (18) 

The integration yields the asymptotic behavior 

D{A)/N (X A-i, (19) 

universal for those threshold distributions for which the assumption of a single 
maximum of {F){x) is valid. 

Note that if the experiment had been stopped before complete breakdown, 
at a force per fiber x less than Xc, the asymptotic behavior would have been 
dominated by an exponential fall-off rather than a power law: 

D{A)/N cx /i-le-[«(^)-i-i"«(^)l^. (20) 

When X is close to Xc the exponent is proportional to {xc — x^A. The 
burst distribution then takes the scaling form 

D{A)ozA-'^G{A''{xc-x)), (21) 

with a Gaussian function G, a power law index /y = | and ^ = \- Thus the 
breakdown process is similar to critical phenomena with a critical point at 
total breakdown [5, 8, 27] . 



2.2 Burst distribution: Nongeneric cases 

What happens when the average strength curve, {F){x), does not have a 
unique maximum? There are two possibilities: (i) it has several parabolic 
maxima, or (ii) it has no parabolic maxima. 

When there are several parabolic maxima, and the absolute maximum does 
not come first (i.e. at the lowest x value), then there will be several avalanche 
series each terminating at a local critical point with an accompanying burst of 
macroscopic size, while the breakdown of the bundle occurs when the absolute 
maximum is reached [9]. The power law asymptotics (19) of the avalanche 
distribution is thereby unaffected, however. 
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The second possibility, that the average strength curve has no paraboHc 
maxima, is more interesting. We present here two model examples of such 
threshold distributions. 

The threshold distribution for model I is, in dimensionless units, 

for a; < 2 

^^""^ ~ \ 1 - (x - for x>2, ^^^> 

while model II corresponds to 
with a positive parameter a. 
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For the two models the corresponding macroscopic bundle strength per 
fiber is, according to Eq. (16), 

(F) _( X for X < 2 

~N~{7§=T^ovx>2 ^^^^ 



for model I, and 



for model II (see Fig. 5). 



(F) ^ / X for a; < 1 , . 

N I x^-" for a; > 1 ^ ' 



To calculate the avalanche distribution we use (16), in both cases with 
the upper limit in the integration. For Model 1 the (F) graph 
has at x = 2+ an extremum, viz. a minimum. At the minimum we have 
a{x) = 2{x-i) ~ ^' the Zi-dependent factor a^e"""^ in the integrand of 
(16) has a maximum for a = 1, we obtain 

DiA)/N (X A-^/'^ (26) 

for large A. Even if the macroscopic load curve does not have a maximum at 
any finite x in this case, the generic power law (19) holds. 
For model II equation (16) gives 

D{A) l-aZ\'^-\ ^„.zi r 

-AT ] ' ] • (27) 

For a = 3/4 as in Fig. 5, or more generally a < 1, the avalanche distribution 
does not follow a power law, but has an exponential cut-off in addition to a 
dependence. 

When a — > 1, the average force (25) approaches a constant for x > 1, and 
the burst distribution (27) approaches a power law 

^ oc A-'/\ (28) 

a result easily verified by simulation on the system with P{x) — 1 — 1/x for 
a; > 1. That a power law with exponent 3/2, different from the generic burst 
distribution (19), appears when 

-i(F)^0, (29) 



will be apparent when we in Sec. 2.5 study what happens at criticality. 
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Let Fk be the force on the bundle when the A:th fiber fails. It is the nonmono- 
tonicities in the sequence Fi, F2, . . . that produce avalanches of size A > 1. Let 
us consider the probability distribution of the force increase AF = Fk+i — Fk 
between two consecutive bursts, the first taking place at a force x = Xk per 
fiber, so that Fk = {N - k + l)x. 

Since AF = [N — k){xk+i — x) — x, it follows that 



AF > -X. 



(30) 



The probability to find the k + 1th threshold in the interval (a;fe+i,a;fe+i + 
dxk+i), for given x = Xk, equals 



(AT 



, ^ Jl-P(a:fc+i)]^-fe-^ , , , 



[1-Pix)] 



(31) 



By use of the connection Xk+i = x+{AF+x)/{N—k) this probabihty density 
for Xk+i is turned into the probability density p{AF) dAF for AF: 

p^AF) = ^'^-^ [1 - P{x +{AF + xyj^N - fe))]^-^-^ ^ ^AF + x 



N-k 



[1-P{x)] 



N^k 
(32) 

which is properly normalized to unity. For large N — k this simplifies to 



piAF) 







l-P(x) 



exp 



(AF+x) 
' l-P(x) 



p{x) 



for AF < -X 
for AF > -X. 



(33) 



P (AF) 




AF 



Fig. 6. The probability distribution p{AF) of the step length in the random walk. 



The values ^1,^2,^3... of the force on the bundle can be considered as 
the positions of a random walker with the probability p{AF) for the length of 
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the next step [8]. It is a random walk of an unusual type: The step length is 
variable, with the steps in the negative direction are limited in size (Fig. 6). 
In general the walk is biased since 

(AF) = I AF p{AF) dAF = ^ " ^^''] ~ ""^^^^ (34) 
J P[x) 

is zero, e.g. unbiased, only at the critical threshold Xc, given by Eq. (8). That 
the random walk is unbiased at criticality is to be expected, of course, since 
the average bundle strength (F) as function of x is stationary here. 
The probability that AF is positive equals 

POO 

PToh{AF > 0) = / p{AF) dAF = exp 

Jo 

That AF is positive implies that the burst has the length A = 1. The result 
(35) is identical to the previously determined probability 0(1, a;), (12), for a 
burst of length 1, when we have not ensured that the burst actually starts 
with the fiber in question and is not part of a larger avalanche. 

In section 2.5 we will see that the random- walk analogy can be used in a 
quantitative way to predict the avalanche distribution power-law exponent at 
criticality. 

2.4 Crossover behavior near criticality 

10° 

10--' 

D(A) 
D(l) 10^ 

10-3 
10" 

1 10 

A 

Fig. 7. The distribution of bursts for throsliold's uniformly distributed in an interval 
{xo,Xc), with xo = and with xo = 0.9xc. The figure is based on 50 000 samples, 
each with N = 10® fibers. 

When all fiber failures are recorded we have seen that the burst distribu- 
tion D{A) follows the asymptotic power law D oc A~^/'^. If we just sample 
bursts that occur near criticality, a different behavior is seen [10, 11] . As 
an illustration we consider the uniform threshold distribution, and compare 



xp{x) 



1 - P{x) 



(35) 
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the complete burst distribution with what one gets when one samples merely 

bursts from breaking fibers in the threshold interval (Q.9xc, Xc)- Fig. 7 shows 
clearly that in the latter case a different power law is seen. 

If we want to study specifically the contribution from failures occurring 
when the bundle is nearly critical, we evaluate the expression (18) for the 
burst distribution over a small interval (a:o, Xc), rather than integrating from 
to Xc- The argument in Sec. 2.1 that the major contribution to the integral 
comes from the critical neighborhood is still valid. We obtain 



N A\ a'{xe) 

with 



1 - e-^/^= 



(36) 



^o = Tr7ZrT27Z-—T2- (37) 



2 

a'{xc)'^{xc - 

By use of Stirling approximation A\ ~ A'^e^^^/2'KA., - a reasonable ap- 
proximation even for small A - the burst distribution (36) may be written 



N 

with a nonzero constant 



= (l - e-^/^») , (38) 



C = {2T,)-^'^p{xc)/a'{xc). (39) 

We see from (38) that there is a crossover at a burst length around A^., so 
that 

^""[a-^I^ for Z\»Zi, (40^ 

The difference between the two power-law exponents is unity, as suggested by 
Sornette's "sweeping of an instability" mechanism [12]. Such a difference in 
avalanche power law exponents has been observed numerically by Zapperi et 
al. in a fuse model [13]. 

We have thus shown the existence of a crossover from the generic asymp- 
totic behavior D cx Z\~^/^ to the power law D oz Zi^"^/^ near criticality, i.e., 
near global breakdown. The crossover is a universal phenomenon, indepen- 
dent of the threshold distribution p{x). In addition we have located where the 
crossover takes place. 

For the uniform distribution Ac = (1 — xo/xc)~'^ /2, so for xq = 0.8 xc, 
we have Ac = 12.5. For the WeibuU distribution P{x) = 1 — exp(— (a: — 1)^°), 
where 1 < a; < oo, we get Xc = 1.72858 and for xq = 1.7, the crossover point 
will be at Ac ^ 14.6. Such crossover is clearly observed (Fig. 8) near the 
expected values A = Ac = 12.5 and A = Ac = 14.6, respectively, for the 
above distributions. 
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A A 

Fig. 8. The distribution of bursts for the uniform threshold distribution (left) with 
xo = O.SOxc and for a WeibuU distribution (right) with xo = 1 (square) and xo = 1-7 
(circle). Both the figures are based on 50000 samples with A'^ — 10^ fibers each. The 
straight lines represent two different power laws, and the arrows locate the crossover 
points Ac ~ 12.5 and Ac ~ 14.6, respectively. 

The simulation results shown in the figures are based on averaging over a 
large number of fiber bundles with moderate N. For applications it is impor- 
tant that crossover signals are seen also in a single sample. We show in Fig. 
9 that equally clear power laws are seen in a single fiber bundle when N is 
large. 

10^ 




1 10 100 

A 

Fig. 9. The distribution of bursts for the uniform threshold distribution for a sin- 
gle fiber bundle with lO'^ fibers. Results with xo = (recording all avalanches), are 
shown as squares, the circles stand for avalanches near the critical point (xo = 0.9xc)- 

An important question in strength considerations of materials is how to 
obtain signatures that can warn of imminent system failure. This is of utter- 
most importance in, e.g., the dimond mining industry where sudden failure of 
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the mine can be very costly in terms of lives. These mines are under continuous 
acoustic surveillance, but at present there arc no tcU-talc acoustic signature 
of imminent catastrophic failure. The same type of question is of course also 
central to earthquake prediction. The crossover seen here in our fiber bundle 
models is such a signature, it signals that catastrophic failure is imminent. 
The same type of crossover phenomenon is also seen in the burst distribution 
of a two-dimensional model of fuses with stochastically distributed current 
trcsholds [10]. This signal docs not hinge on observing rare events, and is 
seen also in a single system (Fig. 9). It has, therefore, a strong potential as 
a useful detection tool. It is interesting that most recently, Kawamura [14] 
has observed a decrease in exponent value of the local magnitude distribu- 
tion of earthquakes as the mainshock is approached (See Fig. 20 of Ref.[14]), 
analysing earthquakes in Japan (from JUNEC catalog). 

Obviously, one cannot count bursts all the way to complete breakdown 
to have a useful detection tool. It suffices to sample bursts in finite intervals 
{xo,Xf), with Xf < Xc- In this case we obtain the avalanche distribution by 
restricting the integration in Eq.(18) to the appropiate intervals. When such 
an interval is in the neighborhood of Xc we obtain 

oc Zi~^/^ ^Q-^i^o-xf f /a _ ^-A(xo-xo)/a^ ^ (42) 

with a = 2/a'{xc)'^. This shows a crossover: 

D{A) r Zi-3/2 for A « a/{xc - xof ^.n. 



N ^\zi-5/2 hva/{xc-xof «A«a/{xc-Xff, 

with a final exponential behavior when A ^ a/{xc — Xf)"^ . 

The 3/2 power law will be seen only when the beginning of the interval, 
Xq, is close enough to the critical value Xc- Observing the 3/2 power law is 
therefore a signal of imminent system failure. 



2.5 Avalanche distribution at criticality 

Precisely at criticality (xo = Xc) the crossover takes place at Ac ~ oo, and 
consequently the ^ = 5/2 power law is no longer present. We will now ar- 
gue, using the random walk representation in section 2.3, that precisely at 
criticality the avalanche distribution follows a power law with exponent 3/2. 

At criticality the distribution (33) of the step lengths in the random walk 
simplifies to 

rATP^ i ior AF < -Xc 

= I a^-^e-^e-^^/- for AF > -Xc ^^^^ 
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A first burst of size A corresponds to a random walk in which the position 

after each of the first A—1 steps is lower than the starting point, but after step 
no. A the position of the walker exceeds the starting point. The probability 
of this equals 



Prob(zi)=r Pc{Si)d6i I ' Pc{S2)d52 I ' \c{S3)dS3... 

J —Xc J —Xc 

-(5i — 52...— <5a — 2 poo 

Pc{5A-i)d5A-i / Pc(^zi)c?^a(45) 

J —5\ —(52 ..-—5^—1 



/"' 

J —Xc 



To sinipHfy the notation wc have introduced (5„ = In Ref. [11] we have 

evaluated the multiple integral (45), with the result 



Prob(/i) 



A^-^e-^ 
A\ 



2tt 



^-3/2 _ 



(46) 



We note in passing that for the standard unbiased random walk with 
constant step length we obtain a different expression for the burst probability, 
but again with a 3/2 power law for large A: 



Prob(Z\) 



1 



2^-1 Z\ V 5^ 



1 



2tt 



^-3/2_ 



(47) 



At completion of the first burst, the force, i.e., the excursion of the random 
walk, is larger than all previous values. Therefore one may use this point as 
a new starting point to find, by the same calculation, the distribution of the 
next burst, etc. Consequently the complete burst distribution is essentially 
proportional to A~^/'^, as expected. The simulation results exhibited in Fig. 
10 are in excellent agreement with these predictions. 



10' c 



D(A) 




Fig. 10. Distributions of the first bursts (squares) and of all bursts (circles) for 
the uniform threshold distribution with xq = Xc- The simulation results are based 
on 10® samples with 80000 fibers each. The crosses stand for the analytic result (46). 
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One of the unusual threshold distributions we studied in Sec. 2.2 corre- 
sponded to an constant average force, (F) /N independent of x. Such a bundle 
is not critical at a single point Xc, but in a whole interval of x. That the burst 
exponent for this model takes the critical value 3/2 is therefore no surprise. 



2.6 Recursive dyneimics 

The relation between the number of ruptured fibers and a given external load 
per fiber, a = F/N , can be viewed as the result of a sequential process. In the 
first step all fibers with thresholds less than xi = a must fail. Then the load 
is redistributed on the surviving fibers, which forces more fibers to burst, etc. 
This starts an iterative process that goes on until equilibrium is reached, or 
all fibers rupture [16, 17, 18]. 

Assume all fibers with thresholds less than Xt break in step number t. The 
expected number of intact fibers is then 

Ut = l- P{xt), (48) 

so that the load per fiber is increased to a/Ut- In step number t+1, therefore, 
all fibers with threshold less than 

must fail. This iteration defines the recursive dynamics [16, 17, 18]. Alterna- 
tively an iteration for the Ut can be set up: 

Ut+i = 1 - P{a/Ut); Uo = 1. (50) 

If the iteration (49) converges to a finite fixed-point x* , 

lim Xt = X*, 

t—>oo 

the fixed-point value must satisfy 
This fixed-point relation, 

a = X* [I - P{x*)], (52) 

is identical to the relation (7) between the average force per fiber, (F) /N, and 
the threshold value x. 

Eq. (52) shows that a necessary condition to have a finite positive fixed- 
point value X* is 

a < (Tc = max {a;[l - P(a;)]} (53) 
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Thus <Tc is the critical value of the external load per fiber, beyond which the 
bundle fails completely. 

As a simple example take the uniform threshold distribution (2) with Xr = 
1, for which the iterations take the form 

xt+i = -7-^— and Ut+i = 1 - (54) 

i- — Xt Ut 

Moreover, cTc = 1/4, and the quadratic fixed-point equation for U* has the 
solution 

U*{a) = i ± (a, - = ± (a, - (55) 

since U*{ac) = 1/2. The iterations are sketched in Fig. 11. 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



u, 

Fig. 11. Graphical representation of the iterations for x (A) and for U (B). The 
graphs are shown for different values of the external stress: a = 0.2 (solid), 0.25 
(dashed), and 0.28 (dotted), respectively. The intersections between the graph and 
the straight 45° line define possible fixed points. 

A fixed point is attractive if \dUt+i/dUt\ is less than 1 at the fixed point 
and repulsive if the derivative exceeds unity. We sec in Fig. 11 that the stable 
fixed point corresponds to the smallest value of x* , and to the largest value 
of U* (the plus sign in (55)): 

U*ia)-U*iac) = {a,-af, with /3 = i. (56) 

Thus U*{a) — U*{ac) behaves like an order parameter, signalling total bundle 
failure when it is negative, partial failure when it is positive. 

Close to a stable fixed point the iterated quantity changes with tiny 
amounts, so one may expand in the difference ct = Ut — U*. To first order 
(54) yields 
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Thus the fixed point is approached monotonously, with exponentially decreas- 
ing steps: 

et (X e-'l\ (58) 

with 

Precisely at the critical point, where U* = \/2 and ac = 1/4, the relaxation 
parameter r is infinite, signalling a non-exponential approach to the fixed 
point. Close to the critical point one easily shows that 

T oc ((Tc — cr)"", with ct=^. (60) 

One may define a breakdown susceptibility x by the change of U*(a) due 
to an infinitesimal increment of the applied stress a, 

^ = -^^ = 5(^c-a)-^ with 7=i. (61) 

The susceptibility diverges as the applied stress a approaches its critical value. 
Such a divergence was noted in previous numerical studies [19, 20]. 

When at criticality the approach to the fixed point is not exponential, 
what is it? Putting Ut = Uc + Cf in the iteration (54) for a = 1/4, it may be 
rewritten as follows 

+ 2, with eo = i, (62) 
with solution e^^ = 2t -|- 2. Thus we have, exactly, 

Ut = - + . (63) 

* 2 2i -I- 2 ^ ' 

For large t this follows a power-law approach to the fixed point, Ut—Uc = ^t~^ , 
with (5 = 1. 

These critical properties are valid for the uniform distribution, and the 
natural question is how general the results are. In Ref. [18] two other threshold 
distributions were investigated, and all critical properties, quantified by the 
indicies a, f3, 7 and S were found to be the same as for the uniform threshold 
distribution. This suggests strongly that the critical behavior is universal, 
which we now prove. 

When an iteration is close to the fixed point, we have for the deviation 

e,+i = -U*=p(-^)-P = • ^P(a/t7*), (64) 

to lowest order in ej. 

This guarantees an exponential relaxation to the fixed point, 
with parameter 
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Critical! ty is determined by the extremum condition (8), which by the relation 
(48) takes the form 

Ul = ap{a/Uc) 

Thus T = oo at criticality. To study the relaxation at criticality we must 
expand (64) to second order in ej since to first order we simply get the useless 
equation ct+i = £*• To second order we obtain 

et+i = e* — Ccf, 
with a positive constant C. This is satisfied by 

Hence in general the dominating critical behavior for the approach to the fixed 
point is a power law with 5 = 1. The values a = /? = 7 = 5 can be shown to 
be consequences of the parabolic maximum of the load curve, (7) or (49), at 
criticality: \xc — cx (cc — cr)2. 

Thus all threshold distributions for which the macroscopic strength func- 
tion has a single parabolic maximum, are in this universality class. 

3 Fiber bundles with local load redistribution 

The assumption that the extra stress caused by a fiber failure is shared equally 
among all surviving fibers is often unrealistic, since fibers in the neighborhood 
of the failed fiber arc expected to take most of the load increase. One can 
envisage many systems for such local stress redistributions. A special case is 
the model with a one-dimensional geometry where the two nearest-neighbor 
fibers take up all extra stress caused by a fiber failure (Sec. 3.1). It is special for 
two reasons: It is an extreme case because the range of the stress redistribution 
is minimal, and, secondly, it is amenable to theoretical analysis. In other 
models, treated in Sec. 3.2 and 3.3, the stress redistribution occurs over a 
larger region. In Sec. 3.3 this comes about by considering a clamp to be an 
elastic medium. 

3.1 Stress alleviation by nearest neighbors 

The simplest model with nearest-neighbor stress redistribution is one-dimensional, 
with the fibers ordered linearly, with or without periodic boundary condi- 
tions. Thus two fibers, one of each side, take up, and divide equally, the extra 
stress caused by a failure. The forcx; on a fiber surrounded by ni broken fibers 
on the left-hand side and rir broken fibers on the right-hand side is then 
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^{l + lni + ^nr)=f{2 + ni + nr), (66) 

where Fioi is the total force on the bundle, and / = F/,o/ /27V, one-half the 
force-per-fiber, is a convenient forcing parameter. The model has been dis- 
cussed previously in a different context[21, 22, 23, 24, 25, 26]. Preliminary 
simulation studies [27, 28] showed convincingly that this local model is not 
in the same universality class as the equal-load-sharing fiber bundles. For the 
uniform threshold distribution and for 1 < Z\ < 10 an effective exponent 
between 4 and 5 was seen, much larger than 5/2. 

Avalanches in this model, and in similar local stress-redistribution models, 
have a character different from bursts in the equal-load-sharing models. In the 
present model the failure of one fiber can by a domino effect set in motion a 
fatal avalanche: If the failing fiber has many previously failed fibers as neigh- 
bors, the load on the fibers on each side is high, and if they burst, the load 
on the new neighbors is even higher, etc., which may produce an unstoppable 
avalanche. 

In Ref.[7] the burst distribution was determined analytically for the uni- 
form threshold distribution, 

P(a:) = |^ forO<x<l 
^ ^ \ 1 for a; > 1 ^ ^ 

In this model there is an upper limit to the size A of an avalanche that 
the bundle can survive. Since the threshold values are at most unity, it follows 
from (66) that if 

^ > - 2, (68) 

then the bundle breaks down completely. Consequently an asymptotic power 
law distribution of the avalanche sizes is not possible. 

In the analytic derivation periodic boundary conditions were used. The 
fairly elaborate procedure was based on a set of recursion relations between 
configurations and events at fixed external force. In Fig. 12 we show the 
resulting burst distribution for a bundle of N = 20000 fibers, compared with 
simulation results. The agreement is extremely satisfactory. 
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Fig. 12. Burst distribution in a local-load-sharing model for a bundle of A'^ = 20000 
fibers. Theoretical results are shown as circles, and simulation results (for 4000000 
samples) are shown as crosses. The straight line shows the power law A~^. 

We see from Fig. 12 that, as expected, the burst distribution does not 
follow a power law for large A, but falls off faster. For A less than 10 the 
burst distribution follows approximately a power law, D{A) oc A~^ with ^ of 
the order of 5. 

In Ref.[7] the maximum load -Fmax the fiber bundle could tolerate was 
estimated to have the following size dependence, 



oc 



N 



(69) 



This is different from the equal-load-sharing model, for which Fmax oc A''. 
A similar logarithmic size-dependence of the bundle strength for local-stress- 
redistribution models has been proposed by other authors [29, 30, 31]. 

The qualitative explanation of the non-extensive result is that for large 
N the probability of finding a weak region somewhere is high. Since, as dis- 
cussed above, a weak region is the seed to complete bundle failure, it may 
be reasonable that the maximum load the bundle can carry does not increase 
proportional to A'', but slower than linear. 



3.2 Intermediate load-sharing models 



It might be interesting to study models that interpolates between global and 
nearest-neighbor load sharing. The main question is whether the burst dis- 
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tribution changes from one behavior to the other in a continuous manner, or 
whether a discontimious change occurs. 

In Ref.[8] was introduced such an intermediate model, with the same one 
dimensional geometry as the nearest-neighbor model of the preceding section. 
When a fiber i fails in this model the elastic constants of the two nearest 
surviving neighbors I and r on both sides are updated as follows 

Ki iA(K; + Kr + Ki) (70) 

Ki ^ (71) 

Kr ^A(«;; + Kr + Ki) . (72) 

For A = 1 this corresponds to the local load-sharing by surviving nearest- 
neighbors (see the preceding section). And with A = the intact neares- 
neighbor fibers to a failing fiber does not take part in the load-sharing. But 
since all the other surviving fibers then share the load equally, this limit- 
ing case must have the same behavior as the equal-load-sharing model. The 
numerics seems to suggest that there is a cross-over value of A separating 
the universality classes of the local-load-sharing and the equal-load-sharing 
regimes. 

A stress redistribution scheme that in a straightforward way interpolates 
between the two extreme models was recently proposed by Pradhan et al. [32] : 
A fraction g of the extra load caused by a fiber failure is shared by the nearest 
neighbors, and the remaining load is shared equally among all intact fibers. 
They show that in a one-dimensional geometry a crossover value gc exists, 
such that for g < gc the bundle belongs to the equal-load-sharing regime, 
while for g > gc the system is like the local-load-sharing model of Sec. 3.1. 
The crossover value was determined to be = 0.79 ± 0.01. 

It would be more realistic to have a stress redistribution whose magnitude 
falls off monotonically with the geometric distance r from the failed fiber. 
Hidalgo et al. [33] introduced such a model, for which the extra stress on a 
fiber followed a power law decay, proportional to r~'^. In the limit 7 — > the 
equal-load-distribution model is recovered, while the limit 7 ^ 00 corresponds 
to the nearest-neighbor model in Sec. 3.1. Again a crossover is observed, at a 
value 7c — 2 for the range parameter. 

In the next section we consider a model with a difii'erent, but similar, 
interaction decaying with the distance from the failed fiber. 

3.3 Elastic medium anchoring 

In this section we generalize the fiber bundle problem to include more realis- 
tically the elastic response of the surfaces to which the fibers are attached. So 
far, these have been assumed to be infinitely stiff for the cqual-load-sharing 
model, or their response has been modeled as very soft, but in a fairly unreal- 
istic way in the local-load-sharing models, see Section 3.2. In [34], a realistic 
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model for the elastic response of the clamps was studied. The model was pre- 
sented as addressing the problem of failure of weldings. In this context, the 
two clamps were seen as elastic media glued together at a common interface. 
Without loss of generality, one of the media was assumed to be infinitely stiff 
whereas the other was soft. 

The two clamps can be pulled apart by controlling (fixing) either the ap- 
plied force or the displacement. The displacement is defined as the change 
in the distance between two points, one in each clamp positioned far from 
the interface. The line connecting these points is perpendicular to the av- 
erage position of the interface. In our case, the pulling is accomplished by 
controlling the displacement. As the displacement is increased slowly, fibers 
— representing the glue — will fail, ripping the two surfaces apart. 

The model consists of two two-dimensional square LxL lattices with peri- 
odic boundary conditions. The lower one represents the hard, stiff surface and 
the upper one the elastic surface. The nodes of the two lattices are matched 
{i.e. there is no relative lateral displacement). The thresholds of the fibers 
are taken from an uncorrclatcd uniform distribution. The spacing between 
the fibers is a in both the x and y directions. The force that each fiber is 
carrying is transferred over an area of size to the soft clamp: As the two 
clamps arc separated by controlling the displacement of the hard clamp, D, 
the forces carried by the fibers increase. As for the fiber bundle models studied 
in the previous sections, when the force carried by a fiber reaches its breaking 
threshold, it breaks irreversibly and the forces redistribute. Hence, the fibers 
are broken one by one until the two clamps are no longer in mechanical con- 
tact. As this process is proceeding, the elastic clamp is of course deforming to 
accomodate the changes in the forces acting on it. 

The equations governing the system are as follows. The force, /», carried 
by the ith fiber is given by 



where k is the spring constant and is the deformation of the elastic clamp 
at site i. All unbroken fibers have k = 1 while a broken fiber has k = 0. The 
quantity (uj — D) is, therefore, the length fiber i is stretched. In addition, a 
force applied at a point on an elastic surface will deform this surface over a 
region whose extent depends on its elastic properties. This is described by the 
coupled system of equations, 



k{ui - D) , 



(73) 




(74) 



where the elastic Green function, Gij is given by [35, 36] 




(75) 
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In this equation, s is the Poisson ratio, e the elastic constant, and the de- 
nominator \i — j\ is the distance between sites i ~ (x,y) and j = {x\y'). 
The indices i and j run over all sites. The integration over the area is 
done to average the force from the fibers over this area. Strictly speaking, the 
Green function applies for a medium occupying the infinite half-space. How- 
ever, with a judicious choice of elastic constants, we may use it for a finite 
medium if its range is small compared to L, the size of the system. 
By combining equations (73) and (74), we obtain 



where we are using matrix-vector notation. I is the L? x identity matrix, 
and G is the Green function represented as an x dense matrix. The 
constant vector D is dimensional. The diagonal matrix K is also x L^. 
Its matrix elements are either 1, for unbroken fibers, or for broken ones. Of 
course. K and G do not commute. 

Once equation (76) is solved for the force /, equation (74) easily yields 
the deformations of the elastic clamp. 

Equation (76) is of the familiar form Kx = b. Since the Green function 
connects all nodes to all other nodes, the x matrix A is dense which 
puts severe limits on the size of the system that may be studied. 

The simulation proceeds as follows: We start with all springs present, each 
with its randomly drawn breakdown threshold. The two clamps are then pulled 
apart, the forces calculated using the Conjugate Gradient (CG) algorithm 
[37, 38], and the fiber which is the nearest to its threshold is broken, i.e. the 
corresponding matrix element it the matrix K is zeroed. Then the new forces 
are calculated, a new fiber broken and so on until all fibers have failed. 

However, there are two problems that render the simulation of large sys- 
tems extremely diflBcult. The first is that since G is a x dense matrix, 
the number of operations per CG iteration scales like L^. Even more serious is 
the fact that as the system evolves and springs are broken, the matrix (I-|-A;G) 
becomes very ill-conditioned. 

To overcome the problematic scaling of the algorithm, we note that 
the Green function is diagonal in Fourier space. Consequently, doing matrix- 
vector multiplications using FFTs the scaling is much improved and goes like 
L^ln(L). Symbolically, this can be expressed as follow: 



where F is the FFT operator and F"-*- its inverse (F~-'-F = 1). Since I and 
K are diagonal, operations involving them are performed in real space. With 
this formidation, the number of operations/iteration in the CG algorithm now 
scales like ln(L). 

To overcome the runaway behavior due to the ill- conditioning we need to 
precondition the matrix [37, 39]. This means that instead of solving equation 
(77), we solve the equivalent problem 



(I + KG)/ = KD 



(76) 



(I + KF-^FG)F~^F/ = KD , 



(77) 



26 



Per C. Hemmer, Alex Hansen, and Srutarshi Pradhan 



Q(I + KF-iFG)F-iF/ = QKD 



(78) 



where we simply have multiphed both sides by the arbitrary, positive definite 
preconditioning matrix Q. Clearly, the ideal choice is Qo = (I + KG)~^ which 
would always solve the problem in one iteration. Since this is not possible in 
general, we look for a form for Q which satisfies the following two conditions: 
(1) It should be as close as possible to Qo, and (2) be fast to calculate. The 
choice of a good Q is further complicated by the fact that as the system 
evolves and fibers are broken, corresponding matrix elements of K are set to 
zero. So, the matrix (I + KG) evolves from the initial form (I + G) to the 
final one I. 

Batrouni ct al. [34] chose the form 



which is the Taylor scries expansion of Qo = (I + KG)"-'^. For best per- 
formance, the number of terms kept in the expansion is left as a parameter 
since it depends on the physical parameters of the system. It is important to 
emphasize the following points: (a) As fibers arc broken, the preconditioning 
matrix evolves with the ill-conditioned matrix and, therefore, remains a good 
approximation to its inverse throughout the breaking process, (b) All matrix 
multiplications involving G arc done using FFTs. (c) The calculation of Q 
can be easily organized so that it scales like riLp' ln(L) where n is the number 
of terms kept in the Taylor expansion, equation (79). The result is a stable 
accelerated algorithm which scales essentially as the volume of the system. 

Fig. 13 (left) shows the force-displacement curve for a system of size 128 x 
128 with elastic constant e = 10. Whether we control the applied force, F, or 
the displacement, D, the system will eventually suffer catastrophic collapse. 
However, this is not so when e = 100 as shown in Fig. 13 (right). In this 
case, only controlling the force will lead to catastrophic failure. In the limit 
when e ^ oo, the model becomes the equal- load-sharing fiber bimdlc model, 
where F = (1 — D)D. In this limit there are no spatial correlations and the 
force instability is due to the the decreasing total elastic constant of the system 
making the force on each surviving bond increase faster than the typical spread 
of threshold values. No such effect exists when controlling displacement D. 
However, when the elastic constant, e, is small, spatial correlations in the 
form of localization, where fibers that are close in space have a tendency to 
fail consequtively, do develop, and these are responsible for the displacement 
instability which is seen in Fig. 13. 



Q = I - (KG) + (KG)(KG) - (KG)(KG)(KG) + ... 



(79) 
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Fig. 13. Force-displacement curve, 128 x 128 systems with e = 10 (left) and e = 100 
(right). 



We now turn to the study of the burst distribution. We show in Fig. 14 
the burst distributions for e = 10 and e = 100. In both cases we find that 
the burst distribution follows a power law with an exponent ^ = 2.6 ±0.1. It 
was argued in Ref. [34] that the value of £, in this case is indeed 5/2 as in the 
global-load-sharing model. 




Burst Size 



10' 10' 
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Fig. 14. Burst distribution for 128 x 128, for e = 10 (left) and e = 100 (right). The 
slope of the straight lines is —2.5. 
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As the failure process proceeds, there is an increasing competition between 
global failure due to stress enhancement and local failure due to local weak- 
ness of material. When the displacement, D, is the control parameter and e is 
sufficiently small (for example e = 10), catastrophic failure eventually occurs 
due to localization. The onset of this localization, i.e. the catastrophic regime, 
occurs when the two mechanisms are equally important. This may be due to 
self organization [40] occuring at this point. In order to test whether this is 
the case, Batrouni et al. [34] measured the size distribution of broken bond 
clusters at the point when D reaches its maximum point on the F — D char- 
acteristics, i.e. the onset of localization and catastrophic failure. The analysis 
was performed using a Hoshen-Kopelman algorithm [41]. The result is shown 
in Fig. 15, for 56 disorder realizations, L = 128 and e = 10. The result is con- 
sistent with a power law distribution with exponent —1.6, and consequently 
with self organization. If this process were in the universality class of percola- 
tion, the exponent would have been —2.05. Hence, we are dealing with a new 
universality class in this system. 




Fig. 15. Area distribution of zones where glue has failed for systems of size 
128 X 128 and clastic constant e = 10. The straight line is a least square fit and 
indicates a power law with exponent —1.6. 
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